% By Martin M. Andreasen, June 2024. 
% For a projection solution, this function computes closed-form expressions 
% for the impulse response functions using the pruning method when using 
% the following definition of the impulse response function .
%
% impulse(x_t+l) = E_t[w_t|eps_t+1=shockSize for the disturbance]
%                 -E_t[w_t]
%
% Note that this definition is the one behind the 
% so-called Generalized Impulse Reponse Function (GIRF)
%
% OUTPUT:
% GIRFy =impulse response function for y = g(x)
%        dimension = ny * length * 2
%        IRFy(:,:,1) = impulse responses at first order
%        IRFy(:,:,2) = impulse responses at second order
% GIRFx =impulse response function for x = h(x)
%        dimension = nx * length * 2
%        IRFx(:,:,1) = impulse responses at first order
%        IRFx(:,:,2) = impulse responses at second order
% NOTE: if xf = [] and xs = [], then we use the default setting and apply
% the mean values of xf and xs.

function [GIRFy,GIRFx] = GIRFPruning_Proj(gx,gxx,h0,hx,hxx,eta,...
    shockSize,xf,order_app,IRFlength)


% The dimension of the model
[ny,nx] = size(gx);
ne      = size(eta,2);
GIRFy   = zeros(ny,IRFlength,2);
GIRFx   = zeros(nx,IRFlength,2);

% We construct the selection matrix S
S       = diag(shockSize ~= 0);

% Defining matrices
Hxx  = reshape(hxx,nx,nx^2);
Gxx  = reshape(gxx,ny,nx^2);   

% Checking on shockSize
[n1,n2] = size(shockSize);
% ensuring that we have a column vector
if n2 > n1
    shockSize = shockSize';
end
if size(shockSize,1) ~= ne
    error('lenght of shockSize must be equal to ne');
end

%% Computing the mean values for xf 
%invI_hx  = inv(eye(nx)-hx);
invI_hx  = (eye(nx)-hx)\eye(nx);
if isempty(xf)
    xf = invI_hx*h0;
end


%% Impulse response functions at first order
GIRFy1st = zeros(ny,IRFlength);
GIRFxf   = zeros(nx,IRFlength);
for k=1:IRFlength
    GIRFxf(:,k)   = hx^(k-1)*eta*shockSize;
    GIRFy1st(:,k) = gx*GIRFxf(:,k);
end
% Saving output
GIRFy(:,:,1) = GIRFy1st;
GIRFx(:,:,1) = GIRFxf;
if order_app == 1
    return
end

%% Impulse response functions at second order
GIRFy2nd = zeros(ny,IRFlength);
GIRFxs   = zeros(nx,IRFlength);
GIRFxfxf = zeros(nx^2,IRFlength);
consQ    = kron(eta*shockSize,eta*shockSize);
hx_hx    = kron(hx,hx);
I_ne     = eye(ne);
vecI_ne  = reshape(I_ne,ne^2,1);
LA       = (kron(eta*(I_ne-S),eta*(I_ne-S)) - kron(eta,eta))*vecI_ne; %Lambda
% Computing GIRFxfxf
for k=1:IRFlength
    GIRFxfxf(:,k) = kron(hx^k*xf+invI_hx*(eye(nx)-hx^k)*h0,hx^(k-1)*eta*shockSize) ...
        + kron(hx^(k-1)*eta*shockSize,hx^k*xf+invI_hx*(eye(nx)-hx^k)*h0) + ...
            hx_hx^(k-1)*(consQ + LA);
end
% Computing GIRFxs and GIRFy2nd
for k=1:IRFlength
    if k > 1
        GIRFxs(:,k) = hx*GIRFxs(:,k-1) + Hxx*GIRFxfxf(:,k-1);
    end
    GIRFy2nd(:,k) = gx*(GIRFxf(:,k)+GIRFxs(:,k)) + Gxx*GIRFxfxf(:,k);
end
% Saving output 
GIRFy(:,:,2) = GIRFy2nd;
GIRFx(:,:,2) = GIRFxf+GIRFxs;
if order_app == 2
    return
end

